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Abstract 

We  report  on  two  recent  advances  in  the  modelling  of  viscoelastic 
polymers:  (i)  a  new  constitutive  model  which  combines  the  virtual 
stick-slip  continuum  “molecular-based”  ideas  of  Johnson  and  Stacer 
with  the  Rouse  bead  chain  ideas;  {ii)  a  two-dimensional  version  of 
a  model  that  accounts  for  stenosis-driven  shear  wave  propagation  in 
biotissue. 


1  Introduction 

The  mathematical  modelling  of  viscoelasticity  (sometimes  also  loosely  re¬ 
ferred  to  as  hysteresis)  in  materials  using  ideas  from  elasticity  has  attracted 
the  attention  of  a  large  number  of  investigators  over  the  past  century.  Among 
signihcant  contributors  (see  the  many  references  in  [15,  16,  17,  18,  19,  28, 
30,  31,  33,  34,  35])  have  been  some  of  the  true  giants  from  the  fields  of 
engineering  and  material  sciences  including  [1]  the  distinguished  mathemati¬ 
cian  to  whom  this  paper  is  dedicated.  One  of  the  most  widely  used  em¬ 
pirical  models  for  viscoelasticity  in  materials  is  the  Boltzmann  convolution 
law  [10,  17,  18,  19,  35],  one  form  of  which  is  given  in  equation  (1)  below; 
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for  summaries  and  further  references,  see  Chapter  2  of  [19]  as  well  as  [10]. 
This  form  of  model,  when  incorporated  into  force  balance  laws,  results  in 
integro-partial  differential  equations  which  are  most  often  phenomenological 
in  nature  as  well  as  being  computationally  challenging.  One  approach  to 
overcome  both  conceptual  and  computational  challenges  is  founded  on  the 
belief  that  hysteresis  is  actually  a  manifestation  of  the  presence  of  multiple 
scales  in  a  physical  or  biological  material  system  that  is  frequently  modelled 
(and  masked)  with  a  phenomenological  representation  such  as  an  hystere¬ 
sis  integral  for  the  macroscopic  stress-strain  constitutive  law.  A  conceptual 
framework  that  is  fully  consistent  with  the  above  viewpoint  is  based  on  inter¬ 
nal  variable  modelling  and  leads  to  an  efficient  computational  alternative  for 
the  corresponding  integro-partial  differential  equation  models.  In  addition, 
it  provides  a  “molecular”  basis  for  the  models  (for  a  comparison  of  mod¬ 
els  of  viscoelastic  damping  via  hysteretic  integrals  versus  internal  variable 
representations,  see  [10]  and  the  references  therein). 

Our  group’s  interest  in  viscoelasticity  in  polymeric  materials  has  been 
motivated  by  projects  in  our  Industrial  Applied  Mathematics  Program  with 
at  least  two  of  our  industrial  partners:  The  Lord  Corporation  and  Meda- 
coustics,  Inc.  The  collaborations  with  polymer  scientists  and  engineers  at 
Lord  involved  the  dynamic  modelling  of  hlled  rubbers  which  experimentally 
exhibit  both  signihcant  hysteresis  and  nonlinearity  in  tensile  and  shear  de¬ 
formations  as  depicted  in  the  sample  stress-strain  curves  in  Figure  1.  The 
efforts  with  engineers  at  Medacoustics  used  some  of  the  viscoelastic  models 
we  have  investigated  in  attempts  to  understand  the  propagation  of  arterial 
stenosis  induced  shear  waves  in  composite  biotissue  in  a  sensor  development 
and  characterization  project. 

In  some  of  our  earlier  efforts  [7,  12,  13],  the  models  for  hysteretic  damping 
in  elastomers  employed  a  phenomenological  Boltzmann-type  constitutive  law 
of  the  form 

/■*  d 

a{t)  =  ge{e{t)) Coeit)  +  J  Y{t  -  s)—gv  {e{s),e{s))  ds,  (1) 

where  e  is  the  inhnitesimal  strain,  Y  is  the  convolution  memory  kernel,  and 
ge  and  g^  are  nonlinear  functions  accounting  for  the  elastic  and  viscoelastic 
responses  of  the  elastomers,  respectively.  As  explained  in  [8,  12],  our  non¬ 
linear  materials  undergoing  large  deformations  required  the  use  of  finite  (as 
opposed  to  infinitesimal)  strain  theories  [31].  However,  since  the  nonlinear¬ 
ity  between  the  stress  and  hnite  strain  is  an  unknown  to  be  estimated  (using 
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Figure  1:  Experimental  stress-strain  curves  for  (1)  unfilled,  (2)  lightly  filled 
and  (3)  highly  hlled  rubber  in  tensile  deformations. 

inverse  problem  algorithms)  and  since  the  hnite  strain  can  be  expressed  in 
terms  of  known  nonlinearities  as  a  function  of  the  inhnitesimal  strain  (at 
least  in  the  problems  we  consider),  one  can  effectively  formulate  the  problem 
as  one  of  estimating  the  unknown  nonlinearity  between  stress  and  inhnites¬ 
imal  strain  (see  [12]).  Hence  one  can  develop  models  for  stress  in  terms  of 
inhnitesimal  strain.  Our  previous  ehorts  as  summarized  in  [8]  have  shown, 
through  comparison  with  experimental  data,  that  the  best  ht  to  hlled  elas¬ 
tomer  data  occurs  when  and  are  cubic,  along  with  F  as  a  distribution 
of  decaying  exponentials.  We  [3,  9]  subsequently  developed  nonlinear  models 
based  on  stick-slip  “molecular”  ideas  of  Johnson  and  Stacer  [23]  and  Doi  and 
Edwards  [16]  which  resulted  in  a  form  for  g^,  gy  and  Y  that  matched  the 
empirical  hndings  reported  in  [8,  12,  13].  These  models  allow  for  multiple 
relaxation  times  present  in  polymer  strands  of  composite  materials  within  a 
virtual  compartmental  model  of  entangled  chemically  cross-linked/physically 
constrained  system  of  long  chain  “molecules” .  While  accounting  for  multiple 
relaxation  parameters,  these  models  do  not  include  physically  or  chemically 
based  parameters  in  representations  of  the  polymer  strands. 

In  the  space  permitted  here,  we  summarize  two  recent  advances:  {i)  a  new 
constitutive  model  [4]  that  has  been  developed  which  combines  the  virtual 
stick-slip  continuum  “molecular-based”  ideas  of  Johnson  and  Stacer  [23]  with 
the  Rouse  bead  chain  ideas  as  described  in  Doi  and  Edwards  [16];  {ii)  a  two 
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dimensional  version  [5,  27]  of  a  model  that  accounts  for  stenosis  driven  shear 
wave  propagation  in  biotissue. 

The  new  molecular-based  constitutive  model,  in  which  polymer  chains  are 
treated  as  Rouse  type  strings  of  interconnected  beads  (a  reasonable  approx¬ 
imation  for  many  materials),  permits  the  incorporation  of  many  important 
physical  parameters  (such  as  temperature,  segment  bond  length,  internal  fric¬ 
tion,  and  segment  density)  in  the  overall  hysteretic  constitutive  relationship. 
Its  form  is  similar  to  that  developed  in  [8,  9]  and  does  have  the  general  form 
(1)  of  Boltzmann  type,  even  though  the  kernel  is  not  of  convolution  type. 

In  the  discussions  of  the  second  part  of  this  paper  we  employ  an  internal 
variable  formulation  of  Boltzmann  type  hysteresis  laws  to  investigate  the 
propagation  of  stenosis  generated  waves  in  biotissue.  It  is  demonstrated  that 
a  viscoelastic  (as  opposed  to  an  elastic)  formulation  is  important  and  that 
waves  generated  in  a  two-dimensional  cylindrical  geometry  with  inner  radius 
partial  occlusions  can  be  readily  modelled  and  simulated. 

2  A  Stick-Slip/Rouse  Hybrid  Model 

We  give  a  brief  outline  of  the  new  constitutive  model;  complete  details  of 
the  derivation  can  be  found  in  the  report  [4].  We  model  a  polymer  mate¬ 
rial  undergoing  directional  deformation  by  assuming  it  is  composed  of  two 
virtual  compartments  as  depicted  in  Figure  2.  One  compartment  consists 
of  a  constraining  tube  which  is  a  macroscopic  compartment  containing  both 
CC  (chemically  cross-linked)  and  PC  (physically  constrained)  molecules.  The 
other  compartment  is  microscopic  in  nature  and  consist  of  those  PC  molecules 
aligned  with  the  direction  of  the  deformation.  These  molecules  will  at  hrst 
“stick”  to  the  constraining  tube  and  be  carried  along  with  its  motion,  but 
will  very  quickly  “slip”  and  begin  to  “relax”  back  to  a  conhguration  of  lower 
strain  energy.  We  compute  the  contributions  of  both  “compartments”  to  the 
overall  stress  of  this  polymer  material  undergoing  deformations. 

To  accomplish  this  we  must  consider  the  contribution  from  the  constrain¬ 
ing  tube  composed  of  both  non-aligned  physically  constrained  molecules  and 
chemically  cross-linked  molecules,  and  that  of  PC  molecules  aligned  in  the 
direction  of  the  deformation  that  are  initially  entangled  with  molecules  of 
the  tube.  These  aligned  molecules  will  in  time  escape  entanglement  and  be¬ 
come  “free”  molecules  and  will  thus  contribute  to  the  overall  stress  in  two 
distinct  phases:  when  entrapped  and  after  “leaking”  free.  Therefore,  there 
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Euh  apped  Portion  of  PC  molecule 


Escaped  Poitioii  of  PC  molecule 


Figure  2:  PC  molecule  entrapped  by  the  surrounding  constraining  tube. 

are  three  contributions  to  the  stress  of  the  system  the  PC  chain  in 

entanglement,  the  portion  of  the  PC  chain  that  has  escaped  entanglement 
and  the  contribution  due  to  the  constraining  tube.  The  constraining  tube  is 
treated  as  elastic  while  we  use  the  Rouse  formulation  [16]  to  treat  the  aligned 
PC  molecules.  We  denote  the  stress  of  the  portion  of  the  polymer  chain  that 
is  constrained  by  the  surrounding  molecules  as  and  the  stress  of  the 

portion  of  the  polymer  chain  that  has  leaked  out  of  the  constraint  tube  as 
Uapit).  The  total  stress  contribution  of  the  entangled  PC  molecules  is  thus 


We  denote  the  stress  of  the  constraining  tube,  assumed  to  be  elastic,  by 
Thus  the  total  polymer  dependent  stress  is  given  by 

C(«)  =  d”’(i)+dT’(«)-  (2) 

Since  we  assume  that  the  entangled  portion  of  the  PC  molecule  behaves 
as  a  free  molecule  immediately  after  a  deformation,  we  will  use  a  particular 
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form  of  the  Rouse  model  [16]  (a  free  molecular  model) 

(X^(t)X^(t)}  =  (1  -  e-^-‘)  (3) 

in  conjunction  with  a  step-strain  process  (similar  to  the  stick-slip  molecular 
formulation  of  Johnson  and  Stacer  [23])  to  arrive  at  an  appropriate  form 
for  the  hysteretic  term  in  a  Boltzmann-type  stress-strain  law.  We 

note  that  Rouse’s  model  is  based  on  the  assumption  that  the  molecule  is 
composed  of  a  continuum  limit  of  a  hnite  number  of  “beads”  connected  by 
some  spring-like  material  (see  Figure  3).  In  (3)  (  is  the  frictional  constant, 
and  Up,  kp  and  ^p  are  constants  depending  on  the  Boltzmann’s  constant  ks, 
the  temperature  T,  the  bond  length  b  and  the  number  N  of  beads  per  chain 
length  (or  chain  length  N  in  the  continuum  limit).  Here  Ro/j(E)  is  the  finger 
strain  given  in  terms  of  the  conhguration  gradient  E.  The  term  W"(t)  in 
Rouse’s  model  represents  the  direction  of  the  time  dependent  Fourier 
coefficient  for  the  bead  (or  at  location  n  along  the  chain  in  the  continuum 
limit)  at  time  t.  To  calculate  the  contribution  of  the  entangled  portion  of 


Figure  3:  Representation  of  vectors  for  a  bead-spring  polymer  molecule. 

the  molecule  to  the  stress  we  subject  the  molecule  to  a  series  of  step-strains 
at  times  0  =  ■  ■  ■  ,tn,  with  the  additional  assumption  that  during  this 

process  some  of  the  entrapped  molecule  will  “leak”  and  escape  entrapment. 
This  “leaked”  portion  of  the  molecule  will  then  be  considered  free,  so  we 
therefore  assume  the  Rouse-like  expression  (3)  describes  its  motion.  Then 
we  let  the  time  between  each  succeeding  step-strain  go  to  zero  in  order  to 
obtain  a  stress-strain  law  that  describes  each  node’s  contribution  to  the  stress 
after  a  step-strain  is  applied  to  the  molecule. 
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If  a  tensile  deformation  is  performed  on  the  elastomer  along  the  axis 
then  the  expression  for  the  constitutive  law  is  given  by  [4] 


S,(«)  =  >y'£Ht)  +  p 


= 


it) 


a. 


(ve)/ 


(elas)  I 


=  (t)  -  (^xx  (t)  +  it)  -  (^xJ  it)  +  it)  -  it) 


&ckftT 


TT^ 


Y.[-2 

p=i  \P  -to 


2X{s)X'is)  + 


X^is) 


ds 


(1  -7(^))' 

p2 


A^(0) 


A(0) 


-2p^ 

e 


+  {  X 


aI-  « 


where  c  is  the  segment  density,  7  is  the  portion  of  the  PC  molecule  still 
constrained  by  the  tube  (7  depends  on  the  tube’s  diameter  a  as  well  as  the 
friction  constant  (  and  bond  length  b),  tr  (which  depends  upon  ^  and  h) 
and  A  =  1  +  e  (with  e  being  the  tensile  strain)  denotes  the  principal  stretch 
which  represents  the  deformed  length  of  a  unit  cube.  The  term  (A^  — 
accounts  for  the  stress  contribution  of  the  elastic  tube  [6]  where  is  the 
Young’s  modulus  of  elasticity.  As  we  mentioned  earlier,  the  integrand  in  (4) 
is  of  the  form 

y(M)4(B„«  =  y(M)4(A7)-^). 

where  Y  is  not  of  convolution  type  because  of  the  factor  7.  Note  that  all  the 
parameters  in  (4)  have  physical  signihcance  attached  to  them. 


2.1  Experimental  Validation 

To  examine  and  validate  the  model,  we  attempted  [4]  to  £t  it  to  several 
sets  of  experimental  data.  One  of  these,  obtained  from  the  experiments 
conducted  by  Huang,  et  al.,  [25],  involved  applying  a  tensile  strain  to  a 
sample  of  articular  cartilage  and  then  finding  the  stress  within  the  cartilage. 
This  data  set  was  used  to  calibrate  and  validate  (4).  The  strain  used  in 
the  experiment  is  a  ramp  strain;  the  strain  increases  at  constant  rate  until 
a  cessation  time  (f^).  The  experimenters  then  levelled  off  the  strain  to  a 
hxed  proportion,  e^ax,  until  the  experiment  terminates  at  time  tf  =  2000 
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seconds.  For  the  experiment,  e^ax  was  taken  to  be  0.05,  with  a  cessation 
time  of  tg  =  400  seconds.  Thus  the  equation  for  the  strain  function  is  given 
by 


e{t) 


-t  t  <ts 

"fy  "fy  n 


(5) 


We  utilized  the  Matlab  tool  Grabit  [20]  to  extract  data  from  Figure  7  in 
[25],  which  is  denoted  here  by  yd-  The  data  extracted  was  used  to  estimate 
the  composite  parameters  9  =  (a,  b,  c,  as  dehned  below  in  (7)  (we  set  T 
to  be  300  K  and  N  to  be  1000)  using  a  least-squares  criterion 


100 


C{d)  =  V  \  -  yd 


i=i 


(6) 


A  Nelder-Mead  (direct  search)  method  (built  into  the  MatLab  function  fmin- 
search)  will  be  used  to  determines  the  optimal  values  for  the  unknown  pa¬ 
rameters.  After  an  appropriate  initial  condition  is  chosen,  the  optimal  value 
returned  by  the  MatLab  program  is 


d 

1 

(M 

3.695  X  10-^ 

b 

b^c 

4.0242  X  10-2 

c 

— 

c 

4.6266 

.  . 

1.9261  X  102 

where  a,  b  and  are  dehned  to  remove  products  of  parameters  in  the  opti¬ 
mization.  If  these  parameter  values  are  used  in  (3),  then  a  very  satisfactory 
(as  conhrmed  by  the  statistical  analysis  in  [4])  £t  of  the  model  to  the  data 
is  obtained  as  depicted  in  Figure  4. 


Stress  versus  Time 


Figure  4:  Stress  calculation  and  comparative  data  for  cartilage  simulation 
with  ramp  strain. 


3  Stenosis-Driven  Shear  Wave  Propagation 
in  Biotissue 

In  this  section,  we  turn  to  recent  results  on  the  viscoelastic  models  for  prop¬ 
agation  of  stenosis-driven  biotissue  waves  mentioned  in  the  Introduction. 
Specihcally  we  report  on  two  dimensional  models  that  employ  an  internal 
variable  approach  to  model  wave  propagation.  To  motivate  this,  we  recall  [2] 
that  coronary  artery  disease  (CAD)  is  caused  by  atherosclerosis,  the  gradual 
accumulation  of  plaque  along  the  walls  of  an  artery.  This  buildup,  known 
as  a  stenosis,  restricts  the  flow  of  blood,  leading  to  a  decrease  in  the  oxygen 
supply  to  the  heart  muscle.  It  is  well  known  that  arterial  stenoses  produce 
sounds  due  to  turbulent  blood  flow  in  partially  occluded  arteries.  In  princi¬ 
ple,  turbulent  normal  wall  forces  exist  at  and  downstream  from  an  arterial 
stenosis,  exerting  pressure  on  the  wall  of  the  artery  which  then  causes  a 
small  displacement  in  the  surrounding  body  tissue.  Our  goal  is  to  model  the 
propagation  of  the  wave  generated  from  the  stenosis  to  the  chest  wall,  and 
ultimately,  to  create  an  inverse  problem  methodology  which  can  be  utilized 
to  determine  the  location  of  an  arterial  stenosis.  Here  we  will  compare  the 
viscoelastic  model  to  an  elastic  one  as  well  as  present  typical  simulations  for 
a  biologically  motivated  example. 

The  simplihed  physical  geometry  which  we  consider,  which  is  based  on 
physical  experiments  at  Medacoustics,  is  a  cylindrical  gel  mold  as  depicted  in 
Figure  5.  The  synthetic  gel  of  which  the  cylindrical  geometry  is  comprised  has 


Figure  5:  Illustration  of  the  physical  geometry  considered  in  our  modelling 
efforts. 
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material  properties  similar  to  those  of  soft  biological  tissues.  A  surgical  tube 
which  mimics  an  artery  with  stenosis  passes  axisymmetrically  through  the 
center  of  the  mold.  A  source  disturbance  (representative  of  the  disturbance 
caused  by  stenosis)  is  generated  within  the  tube.  Shear  waves  propagate 
through  the  gel,  and  acceleration  is  measured  (by  sensors)  at  the  outer  surface 
of  the  gel.  The  following  system  of  equations  were  derived  from  hrst  principles 
in  [27]  for  this  physical  geometry: 


P 
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d'^ui 

dt"^ 

d‘^U2 

dt"^ 

dex, 

dt 

dt 

dt 

(9e22 
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^(eA  +  €“)  +  +  -(4‘  ~  V 
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+  C'Afc^  +  S'22^j 
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—  +  C'itl.TW  (  “S'. 


d 

Mk 
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12 


i^k 


dt 


;(e) 

-*22 
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(8) 

(9) 

(10) 

(11) 

(12) 

(13) 


where  Ui  and  U2  are  displacement  in  the  radial  and  tangential  directions, 
respectively,  e\^  and  represent  the  K  internal  variables,  z/j  and  Ci  are 
material  parameters.  The  represent  the  elastic  response,  given  by  the 
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—  2a^Ei2  +  2ca^Ei2  exp(aii744  +  027722  +  ®3-^i2 
+037721  +  2a/^EiiE22) 

=  2(q;477ii  +  0:27722)  +  2c(o477ii  +  O27722)  exp(oi77 

+027722  +  0377i2  +  a^E^i  +  2o477ii7722). 

Several  different  variations  of  the  elastic  response  were  considered  in  [27].  For 
the  hrst  elastic  response  model  (ERl),  it  was  assumed  that  the  off  diagonal 
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strain  terms,  E12,  were  equal  to  zero,  and  that  the  small  stress  terms  were 
negligible  (0;^  =  0).  The  second  elastic  response  model  (ER2)  was  formulated 
under  the  assumption  that  the  small  stress  terms  were  negligible,  while  the 
third  elastic  response  model  (ER3)  retained  the  small  stress  terms  and  as¬ 
sumed  that  the  off  diagonal  strain  elements  were  zero.  Statistical  signihcance 
tests  were  employed  to  determine  which  versions  of  the  model  produced  best 
hts  to  radially  symmetric  (and  hence  1-dimensional)  data  presented  in  [2], 
From  this  it  was  determined  that  the  hrst  elastic  response  model  employing 
two  (i.e.,  K  =  2)  internal  variables  (denoted  by  (ER1-2ISV))  provided  one  of 
the  best  hts.  Thus,  the  results  presented  here  were  all  generated  using  that 
model. 

Equations  (8)  and  (9)  are  variations  of  the  equations  of  motion.  They  are 
derived  from  the  physical  properties  of  conservation  of  mass  and  conservation 
of  momentum.  The  remaining  equations  (10)-(13)  represent  the  dynamics  of 
the  internal  variable.  They  are  based  on  the  internal  variable  approach 
employed  in  [2]  as  an  approximation  of  and  alternative  to  Fung’s  quasi-linear 
viscoelastic  theory. 

In  [19],  Fung  develops  and  presents  the  quasi-linear  viscoelastic  constitu¬ 
tive  equation 

s^,{t)  =  £  (14) 

where  Sij  is  the  Kirchoff  stress  tensor,  E  is  the  Green’s  strain  tensor,  Gijki 
is  a  reduced  relaxation  function,  and  is  the  “elastic”  stress  tensor. 

Since  its  introduction,  this  quasi-linear  viscoelastic  (QLV)  theory  has  been 
applied  successfully  in  stress-strain  experiments  to  several  types  of  biological 
tissue.  A  beneht  to  using  (14)  as  a  constitutive  equation  is  that,  unlike  sim¬ 
pler  models  for  viscoelasticity,  it  allows  for  the  consideration  of  a  continuous 
spectrum  (e.g.,  see  the  discussions  in  [19])  of  relaxation  times  and  frequen¬ 
cies  (this  is  also  true  of  the  probabilistic-based  internal  variable  approach  we 
developed  in  [11]).  While  Fung’s  theory  has  been  successfully  employed  for 
htting  hysteretic  stress-strain  curves,  we  are  interested  in  using  it  in  a  full 
dynamical  model.  Unfortunately,  the  QLV,  as  presented  by  Fung,  leads  to 
exceedingly  difficult  computations  within  full  dynamical  partial  differential 
equations,  especially  in  inverse  problems.  This  motivated  the  development  of 
the  internal  variable  approach  described  in  [2,  11,  27]  (which  permits  discrete 
approximation  to  a  continuum)  in  attempts  to  approximate  well  the  corre¬ 
sponding  dynamic  responses  even  in  cases  where  the  stress-strain  curves  alone 
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do  not  produce  adequate  approximations  -  see  [19]  for  discussions. 

Our  alternative  to  Fung’s  kernel  is  a  parameter  dependent  kernel  with  a 
continuous  distribution  of  parameters  and  internal  variables.  In  the  case  of 
a  hnite  combination  of  Dirac  S  distributions  (employed  in  this  presentation), 
we  obtain  a  hnite  summation  of  exponential  functions  as  our  approximation 
kernel.  This  method  can  be  extended  to  allow  for  consideration  of  a  con¬ 
tinuous  spectrum  of  relaxation  times  and  frequencies  by  utilizing  absolutely 
continuous  parameter  distributions  in  place  of  the  6  distributions. 

To  numerically  integrate  our  system,  we  hrst  used  a  series  of  substitutions 
so  that  we  could  rewrite  it  as  a  system  of  hrst  order  partial  diherential  equa¬ 
tions  for  a  vector  U  containing  velocities  Vi  and  strains  tCj.  We  utilized  the 
MacCormack  hnite  diherence  scheme  (see  [33])  to  numerically  integrate  the 
system.  This  scheme  is  a  two-step  method  with  second  order  accuracy  which 
uses  forward  diherences  on  the  hrst  step,  followed  by  a  step  of  backward 
diherences. 

To  update  boundary  terms,  we  used  direction  cosines  (see,  e.g.,  [24]).  If 
the  characteristic  curve  points  into  the  domain,  then  a  physical  boundary 
condition  is  prescribed  for  the  boundary.  Otherwise,  the  boundary  values 
are  updated  utilizing  the  second  order  discretization 

un+i  =  At[yh7,  +  BUg  +  Qr  +  ^At^[AUr  +  BUg  +  Q]l  (15) 

As  there  is  only  one  characteristic  pointing  into  the  domain  at  both  the  inner 
and  outer  radius,  we  prescribed  the  following  boundary  conditions.  On  the 
inner  radius,  wi{Ri,9,t)  =  fit),  where  fit)  is  an  approximate  impulse  input 
function  similar  to  that  used  in  [2].  On  the  outer  radius,  we  prescribe  wi  =  0. 

Although  our  models  (and  those  of  Fung)  have  been  derived  under  the 
widely  accepted  assumption  that  biotissue  has  viscoelastic  properties,  there 
are  other  groups  [21,  22,  29]  conducting  research  in  this  held  that  utilize 
non-dissipative  constitutive  equations  related  to  elastic  media  for  their  mod¬ 
elling  efforts.  This  raises  the  question  as  to  whether  the  more  complicated 
viscoelastic  models,  as  proposed  by  Fung,  our  group  and  others,  really  are 
required.  We  therefore  examined  the  difference  between  a  model  employing 
an  elastic  assumption  and  one  which  assumes  viscoelasticity  as  formulated 
here.  The  Maxwell,  Voigt,  and  Kelvin  models  for  viscoelasticity  are  rather 
easily  reduced  to  an  elastic  model  by  setting  the  viscous  portion  to  be  zero. 
In  a  similar  manner  we  can  easily  reduce  our  model  to  a  purely  elastic  model, 
and  thus  the  models  are  readily  compared. 
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By  considering  a  summation  of  exponential  functions  to  approximate  the 
Fung  kernel,  we  defined 

A(f)  =  fxit)  = 

As  detailed  in  [27],  if  we  allow  the  relaxation  parameters  ux  and  to  tend 
towards  zero,  then  our  model  reduces  to  a  purely  elastic  model.  To  examine 
the  effect  of  the  u  parameters  on  our  system,  we  ran  simulations  with  values 
of  and  Uk  =  0,  where  represents  the  optimized  parameter 

from  [27]  using  the  radially  symmetric  data. 

All  simulations  and  boundary  conditions  were  treated  as  explained  pre¬ 
viously.  Figure  6  depicts  time  snapshot  results  for  vi  for  different  values  of 
V  at  two  different  time  steps. 


Figure  6:  Results  for  vi  for  the  varying  values  of  the  relaxation  parameter, 
V  at  t=0.0032  (left)  and  t=0.0504  (right). 


We  observe  from  Figure  6  that  as  time  increases,  the  solutions  produced 
using  larger  values  for  the  relaxation  parameter  exhibit  more  dissipation.  By 
time  t  =  0.0504  seconds  (depicted  in  Figure  6  (right)),  it  is  evident  that 
the  only  wave  which  has  shown  no  dissipation  is  that  which  was  generated 
using  the  purely  elastic  model  (i.e.,  z/  =  0).  This  is  what  one  would  expect 
from  a  wave  travelling  in  a  purely  elastic  medium.  The  wave  should  simply 
reflect  from  one  boundary  to  the  other,  travelling  back  and  forth  with  no 
dissipation  (which  we  note  does  NOT  agree  with  behavior  in  any  biotissue 
material).  Our  results  illustrate  that  setting  z/  =  0  effectively  reduces  our 
model  from  viscoelastic  to  elastic  in  nature. 
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While  our  results  demonstrate  that  our  model  can  be  easily  reduced  to 
an  elastic  model,  they  also  indicate  that  the  results  produced  by  an  elastic 
model  vary  signihcantly  from  the  results  generated  by  a  model  incorporating 
viscoelasticity.  Since  we  know  dissipation  occurs  in  biotissue,  this  supports 
the  argument  that  the  viscoelastic  model  clearly  is  a  more  realistic  model. 
The  difference  between  the  elastic  and  viscoelastic  simulations  motivate  our 
efforts  in  using  a  viscoelastic  approach  to  modelling  wave  propagation  in 
biotissue. 

Having  numerically  motivated  the  use  of  a  viscoelastic  model,  we  next 
present  simulation  results  generated  by  the  model  for  a  geometry  modihed 
to  represent  a  stenosis.  We  assumed  that  a  buildup  of  plaque  (cholesterol, 
calcium,  and  platelets)  has  formed  along  the  wall  of  the  inner  radius  of  our 
geometry  (see  Figure  7).  We  assumed  that  this  buildup  is  completely  rigid 
and  impermeable,  thus  any  impulse  along  the  inner  radius  will  have  no  direct 
effect  on  the  region  occluded  by  the  buildup. 


Figure  7:  Geometry  with  occlusion  along  inner  radius. 


In  order  to  numerically  solve  our  system  under  these  conditions,  we  shall 
dehne  the  impulse  along  the  inner  radius  as  follows:  Wi(t;  Ri,9)  =  fi(t)  * 
s{9),  where  fi(t)  is  the  impulse  function  used  above  and  s{9)  is  dehned 
as  the  characteristic  function  for  [0,6^1)  U  (6^2, 27r].  For  the  computational 
simulations  presented  here,  we  assumed  that  the  occlusion  occupied  the  area 
between  9i  =  ^  and  6*2  =  ^.  Computations  were  again  carried  out  using 
the  MacCormack  scheme,  and  boundary  values  were  again  updated  using 
direction  cosines.  However,  the  input  function,  Wi(t;  Ri,9),  as  dehned  above. 
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now  has  discontinuities  at  the  edges  of  the  occlusion  (^i  and  6*2).  These 
discontinuities  present  numerical  difficulties  with  the  integration  scheme,  and 
must  be  treated  with  certain  precautions  [24,  26].  Although  the  method 
of  direction  cosines  permits  one  to  prescribe  only  one  boundary  condition 
(corresponding  to  the  one  eigendirection  that  points  into  the  domain)  on 
the  inner  radius,  the  presence  of  discontinuities  in  the  input  permits  the 
assignment  of  additional  conditions.  To  ensure  that  the  problem  is  well 
defined,  we  must  dehne  the  velocity  across  the  discontinuity.  Hirsch  [24] 
argues  that  in  this  situation,  one  should  dehne  the  tangential  velocity  to  be 
zero  across  the  discontinuity.  This  means  that  we  may  prescribe  V2{t]  0,  61  ± 
e)  =  0  and  V2it;  0,  6*2 ±e)  =  0.  We  ran  our  simulations  using  the  discontinuous 
inner  radius  boundary  conditions,  as  well  as  the  new  boundary  conditions; 
typical  results  are  displayed  in  Figures  8  and  9. 


Figure  8:  Results  for  vi  (left)  and  V2  (right)  for  the  geometry  with  an  occluded 
inner  radius,  t=0.0039. 

By  examining  Figures  8  and  9,  we  observe  that  as  time  increases,  the 
radial  velocity  (depicted  by  the  graphs  on  the  left  of  each  hgure)  propagates 
outward  uniformly  throughout  the  geometry,  with  the  exception  of  the  area 
affected  by  the  occlusion.  The  tangential  velocity  is  also  propagating  out¬ 
ward,  but  this  only  occurs  around  the  interface  between  the  occluded  and 
non-occluded  portion  of  the  inner  radius. 

The  results  summarized  here  and  further  detailed  in  [5,  27]  illustrate  an 
ability  to  simulate  stenosis-driven  wave  propagation  through  homogeneous 
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Figure  9:  Results  for  vi  (left)  and  V2  (right)  for  the  geometry  with  an  occluded 
inner  radius,  t=0.0162. 


biotissue.  Of  course,  the  composition  of  the  human  chest  cavity  is  much 
more  complex  than  the  geometry  considered  in  the  model  presented  here. 
The  inclusion  of  heterogeneities  due  to  different  tissues  (lung,  bone,  etc.) 
present  in  the  chest  cavity  is  necessary  for  improvement  of  our  modelling 
efforts.  (Initial  results  in  this  direction  are  presented  in  [5,  27].)  Nevertheless, 
this  model  and  its  simulation  capabilities  provide  a  solid  foundation  for  the 
inclusion  of  heterogeneities,  as  well  as  a  good  starting  point  for  examining 
inverse  problems  formulated  to  determine  the  location  and  extent  of  a  stenosis 
within  a  given  two  or  three  dimensional  geometry. 
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